FLOW FIELD ANALYSIS OF 
TWO DIMENSIONAL SUPERSONIC INTAKE 


A Thesis Submitted 

In Partial Fulfilment of the Requirements 
for the Degree of 

MASTER OF TECHNOLOGY 


by 

VINAYAK SISTA 


to the 

DEPARTMENT OF AEROSPACE ENGINEERING 

INDIAN INSTITUTE OF TECHNOLOGY, KANPUR 

MAY, 1989 



‘3 OC 11989 


CEMr;--L library 

'■•::? ?>o-' 

! V . KANPUH 


icc.No. 

A'-' 

A£ “ f\q - s ^ 3 _ 

FLC7 



CERTIFICATE 


This is to certify that the thesis work entitled, 
” Flow Field Analysis of a Two Dimensional Supersonic Intake” 
has been carried out by Sista. Venkata. Vinayaka.Vishwanatham 
under my supervision and that it has not been submitted 
elsewhere for a degree. 


c 

Dr. O.P. SHARMA 

Professor 

Department of Aeronautical Engg. 
Indian Institute of Technology 


Kanpur . 



CONTENTS 


NOMENCLATURE Page 

CHAPTER 1 : INTRODUCTION ^ 

CHAPTER 2 : FORMULATION OF THE PROBLEM 9 

2.1 : MATHEMATICAL FORMULATION 11 

CHAPTER 3 : METHOD OF SOLUTION 22 

3.1 : TEST CASE 25 

CHAPTER 4 : DISSIPATION 32 

CHAPTER 5 : RESULTS, DISCUSSIONS AND CONCLUSION 38 


LIST OF REFERENCES 


52 



NCaiENCLATORE 


X = Ratio of specific heats = 1.4 
u = Velocity component in x direction 
V = Velocity component in y direction 
p = Pressure 
T = Temperature 

R = Charecteristic gas constant 
U,F,G = Fluxes in conservation form 
ho = Stagnation or total enthalpy 
e = internal energy 
A = Area of the C 

Hab» Hgc, Hcd, Hda = Are the fluxes through the sides 
XA» xg = X coordinates at A and B 
yA» YB - y coordinates at A and B 
w = Relaxation factor 
0 = Wedge angle 
P = Shock angle 

M = Mach numbers ’ 

X = Step length in x direction 
y = Step length in y direction 
t = Time step 
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Chap'ter 1 


GENERAL INTRODUCTION 

Uith the advent o£ auperaonic aircr-a£t , the design o£ the 
intake has become critical as it has a marked e££ect on the 
engine e££iciency and external drag. In the early stages, the 
performance was primarily predicted experimentally, that is, 
from wind tunnel testing. But now the availability of high 
speed, high storage capacity computers, and algorithmic 
sophistication has allowed the application of advanced 
analytical methods to several practical design problems 
including the inlet design and performance analysis. 

Modern computational methods allow direct numerical 
solution of inlet flow field in an efficient manner. A variety 
of integrated performance parameters can be obtained from such 
a solution. Thus different configurations can be developed 
with less reliance on wind tunnel testing, especially in the 
preliminary design stages. However, computational fluid 
dynamics is still a growing field and for too complex flows, 
exact solution is not possible even with present algorithms and 
computing machines. In such cases, exact solutions are obtained 
for simpler configurations to establish correctness of the 
procedure and then modified to handle complicated cases. 

The primary function of an inlet is to supply 
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continuous air flow to the engine at a desired rate with 
highest possible recovery of total pressure, the lowest 
possible external drag and miniinuin weight. Further it should 
also provide uniform flow distribution at the entry of the 
compressor or combustion chamber. Usually the entry Ilach number 
is around 0.4 for subsonic combustion. For a uniform flow 
field and design IHach number, an intake will have maximum 
possible pressure recovery, but for a wide range of mach 
numbers and angles of attack, the performance is far from 
optimum. During the deceleration of airflow, the static 
pressure is increased. Since the process of deceleration can 
be achieved both inside as well as outside the inlet, it is 
convenient to discuss each separately. 

INTERNAL COMPRESSI TYPE 

Here the compression is achieved through a series of 
internal oblique shocks followed by normal shock positioned 
downstream of the throat of the intake. Finally the subsonic 
flow downstream of the constant area throat is further slowed 
down by moving through a divergent section. 

EXTERNAL COMPRESSICMJ TYPE 

The compression in this type of inlet is achieved by 
either one or a series of oblique shocks followed by a normal 
shock. The external compression inlet that achieves compression 
through single normal shock is called Pitot inlet or Normal 
shock inlet. This type of inlet is simple, short and light in 
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weight. Total pressure recovery of this inlet corresponds to 
total pressure ratio across the normal shock. However for 
numbers greater than 1.6^ The total pressure recovery is too 
low and a more efficient design should be used. 

The external compression inlet with one or more oblique 
shocks has its inlet throat at or very near the cowl lip. 
Increasing the number of oblique shocks increases the total 
pressure recovery for any free stream Mach number. The 
compression is achieved by a ramp. This type of inlet can be 
used upto Mach number 2.5 with acceptable pressure ratio and 
cowl drag. 

MIXED COMPRESSION TYPE 

This type of inlet achieves compression through the 
external oblique shocks, internal reflected oblique shocks, and 
the terminal normal shock. Ideal location of normal shock is 
just downstream of inlet throat to minimize total pressure 
losses and to maintain a stable operating location of this 
shock. It can be used above Mach 2.5: It is more complex, 
heavier and costlier than external compression types. 

The supersonic inlets can also be classified as two 
dimensional (rectangular )and axisymmetric (circular). 
Axisymmetric inlets have a alight advantage over two 
dimensional inlets with respect to weight and total pressure 
ratio. However, the two dimensional inlets have an advantage in 
design simplicity and in providing larger variation in inlet 
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the flow field for an arbitrary geometry. The geometry selected 
is as shown in fig(l.l). The upper wall is kept parallel to the 
free stream, whereas at the lower wall the flow undergoes a 10* 
compression at 2 cm. and a 10* expansion at 4 cm. Total length 
of the intake is 10 cm. and height is 2 cm. In this inlet we 
study the variations in pressure, density, and velocity across 
the shocks. From this we know the maximum possible pressure 
recovery of the inlet. It is a complex problem with shock 
expansion layer interaction and shock boundary layer 
interaction etc. So, for sake of simplicity, certain 
assumptions are made. These are described in the next chapter. 
For this problem the solutions are obtained for inflow 
condition as given, Mach number=5, P=1.013, and T=298K. After 
calculations the results are plotted and compared. 

1.2. LITERATURE REVIEW 

In the area of inviscid compressible flows there are a variety 
of practical problems. These include rocket nozzle flows, 
aircraft and missile inlet flows, etc. The methods that have 
evolved and are used currently for the solution of compressible 
flows fall roughly into two categories. The first is the 
unsteady method which is used to integrate the flow eqns. with 
respect to time for both steady and unsteady flows. This 
approach has been useful in overcoming difficulties associated 
with shock waves and transonic flows. The second approach is 
the marching techniques for integration of steady eqns. in 



supersonic flows. 


The use of unsteady methods for supersonic free streams 
is primarily for the computation of unsteady or steady flows 
with embedded subsonic region. Bohachevsky and Rubin (1966) 
applied unsteady finite difference methods to Moo > 1 flows. 
These authors used the highly damped method of Lax (1954), to 
solve inviscid Euler eqns. for flow past a blunt body. Next 
came Horrety and colleagues using schemes suggested by Lax and 
Uendroff (1960 , 1964 ) and, successively by the Maccormack 
method(1969) . The solutions of Lax-Uendrof f , Maccormack and 
Rusnow methods worked satisfactorily for flow past bodies with 
fairly regular shapes such as wedges, spheres and cones, but 
have a tendency to fail for complicated geometries. 

For smooth flows one can use the Maccormack scheme 
(1981) to integrate the inviscid flow eqn. and obtain good 
results. For steady flows one can replace the unsteady energy 
eqn. by steady constant-total- enthalpy eqn. and still obtain 
good results. This was tested for compressible inviscid flows 
for different geometries and fairly good results were obtained. 

Morrety (1978, 1979) set up a method that examined 
characteristic directions and selected difference formulas to 
approximate the flow eqns. depending on behavior of the 
characteristics. His procedure was to write the inviscid flow 
eqns. in the form which in characteristic methods is called 
compatibility eqns. If these are integrated in characteristic 
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direction they reduce to ordinary differential eqns., which are 
much simpler to deal with. He tested this scheme for flow past 
blunt bodies and aircraft configurations successfully. The 
method of characteristics gives a better solution than any 
other method but the drawback is the tedious procedure involved 
in computing. 

IMPLICIT METHODS 

The limitations on the time step induced by stability 
considerations associated with explicit schemes is often too 
restrictive in applications, and consequently implicit schemes 
are used more and more to overcome this limitation. 
Ritchermeyer and Morten (1967) proposed this scheme for the one 
dimensional case. This method has restricted use due to 
programming complexities and computing time and stability 
restrictions. A jay Kumar (1981) developed a new method 
combining explicit and implicit methods and tested it on a 
supersonic intake model. Results were fairly accurate and there 
was saving in computer time. 
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Chapt,er 2 


Forimila'Llon of 'the problem 

The basic assumptions in solving the flow field of a 
supersonic inlet are as follows: 

1. Intake is two dimensional. 

2. Flow is steady, inviscid and adiabatic. 

3. No body forces 

4. Turbulence effects neglected. 

5. Ratios of specific heats ir') and entropy are assumed constant. 

Usually axisymmetric intakes are used due to higher 
pressure recovery but due to design simplicity a two 
dimensional model has been selected. In reality the flow is 
viscous, hence there is formation of boundary layer, heat 
transfer effects near wall and many other undesirable effects. 
In supersonic flows shock formation will be present and in the 
presence of boundary layer, there will be a shock wave boundary 
layer interaction which may lead to separation. In the 
supersonic intakes the temperature also varies rapidly after 
the shock, which leads to variation in Cy ) and rise in entropy. 
The analysis with the above mentioned things is too 
complicated, so the present analysis is restricted to steady 


inviscid flow. 
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FUNDAMENTAL DIFFICULTIES 

Basically there are four dependent variables in any 
ideal fias, two dimensional flow problem: Two velocity 
components and two thermodynamic properties. In the case of an 
incompressible flow, the energy equation is not required to 
solve the momentum and continuity equations, removing one 
thermodynamic property, namely temperature. The pressure may be 
removed by cross differentiation, and vorticity may be 
introduced. The two velocity components are removed by 
introduction of stream function, leaving the two unknowns 
vorticity and stream function to be solved for by two 
equations. In case of compressible flow the energy equation is 
required to solve the others, and stream function is not 
defined for unsteady flow. So we must work with four coupled 
partial differential equations. 

The inviscid equations for supersonic flow are 
hyperbolic. The most formidable numerical problem unique to 
supersonic flow is existence of shock waves. In high 
Reynold’s number limit, shock waves are discontinuities in the 
flow solution. Most of the efforts in this area have been 
directed towards smearing (capturing) out these 
discontinuities . 

The numerical computation is complicated, because the 
regular solution may break down due to nonlinearity of the 
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governing equations. The breakdown is characterized by shock 
wave formation. In the last decade two numerical techniques 
have emerged for the analysis of these flows. One is the shock 
capturing technique. It tries to remove explicit computation of 
discontinuities by generalizing the concept of a solution of 
an Euler equation to include weak solutions (i.e. 
discontinuities). As this scheme requires no special treatment 
to deal with discontinuities, it has become extremely popular 
in computing. However despite its present popularity it is a 
poor interpretation of the physical phenomenon and an 
uneconomical way of computing. 

The other technique is shock fitting. It makes special 
provisions for explicitly computing discontinuities. In essence 
it locates the discontinuities and treats them as boundaries 
between regions where a regular solution is valid. 


2.b MATHEMATICAL FORMULATKW 

The governing equations for an inviscid, two dimensional 
incompressible flow are as follows: 

1. Continuity equation: 


££ + _ 
^t ^x ^y 


2 . 


Momentum eqn. in x-direction 
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du . du . dp 


2.2 


3. 


Momentum eqn. in y direction 


dv . 4 , dv ^ dv dp 

"yt ■* ♦ '=''55 ‘ -55 


2.3 


Energy eqn. 


d . V d V d V 

3t< ® - 2 > * "“55 "''55 < “ 2 ^ 

= -fe < ^ - I 5 < '“ > 


5. Equation of state 

P = PRT ..2.5 

CONSERVATION FORM: 

The conservation form of inviscid compressible flow 
eqns.was given by Courant and Friedrichs (1948) but Lax (1954) 
was the first to use this form to achieve a conservative finite 
difference scheme. 

Uhen the usual differential equations are rearranged so 
that the conservation variables p, pu, pv, e (to be defined) 
are the independent variables, then the use of conservative 
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finite difference methods aasures identical conservation of 
mass, momentum and energy. The Rankine Hugnoit relations (2) 
for a normal shock wave are based on gross conservation and are 
independent of details within the shock structure. The 
conservation equations satisfy the Rankine Hugnoit relations 
and therefore produce correct jump across the shock. Longley 
(1960) tested four separate differencing methods and found that 
all produce correct shock speeds because conservation equations 
were used. Many subsequent calculations showed that 
conservation form is more accurate for calculations across the 
shock wave. 

The conservation form can be obtained by manipulating 

the equations (2.1 2.5). The final form of the equation is as 

follows : 


^ £F ^ £G 
dt dx ^y 


O 


.. 2.6 


Uhere U, F, G are defined as 
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pu 

pu*+ p 

F= 

puv 
puh 

o 

pv 
puv 

G ' I 

f,v * p 

pvh 

Thus we have Euler’s equation (2.6) in conservation form. 

Ref. (5 ) 

Assufiiing that enthalpy remains constant i.e. 

h = C ^ 5 - + — fu^+v*) 27 

o y-1 p 2 ^ ' 

Thus the energy equation is eliminated from the solution and 
pressure can be directly determined from the above eqn. 

Now by knowing the initial and final conditions for the 
flow eqn. (2.6) can be integrated to provide the inviscid 
solution at a later time level. Since steady state is a special 
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case of unsteady flow, the steady flow solution can also be 
obtained by solving the unsteady Euler eqn. 

The eqns . may be solved in either Finite difference form or 
Finite volume form e.g. (Veuillot (1976), Gopalakrishna and 
Bozzala (1971)). It is usual to transform the computational 
domain into uniform rectangular grid and to express derivatives 
of flow variables in terms of values at the nodes of the grid. 
Specialized numerical techniques (Haccormack or Lax-Uendroff 
schemes are needed to ensure the stability of the integration 
of the equation through time until steady state is reached. 

In finite volume method (HcDonald (1971), Denton (1975)) 
the equations are regarded as equations for conservation of 
mass, energy and momentum applied to a set of interlocking 
control volumes formed by a grid in the physical plane. Thus in 
this way the conservation of all the three is ensured. This is 
not true in the differential approach. In the present analysis 
the Finite volume formulation is used. 

FINITE VOLUME FORMULATION 

The Finite volume cell method is based upon an integral 
form of equation to be solved. The computational region is 
divided into elementary quadrilateral volumes within which 
integration is carried out. Such a procedure allows one to deal 
with complicated geometry without considering the equation 
written in curvilinear coordinates. (This also preserves 
the property of conservation). Only the coordinates of the 
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corners of the cell are necessary. The technique is desicribed 
by a model equation. 


^t Sx Sy 


0 


The integral form of the above equation is 


&t 


dA 


/ V.H dA 


= 0 


Assuming 


dU 


constant within the cell, 


du 

dt 


A 


dA 


+ 


/ H.ndS 
S 


0 



+ / C F.l + 

S * 


G.l >dS 

y 


0 


. .2.7 


Uhere A = Area of the cell 

S=surface integral 

n= unit normal to the boundary 

If cartesian coordinates are used in the evaluation of boundary 
integral then we have, 

H.ndS = Fdy - Gdx 

Assuming a curvilinear mesh Fig (2.1), the elementary cell 
is quadrilateral ABCD. The eqn. 2.7 reduces to 


= ( H 


AB 


H + 

BC 


H 


CD 


H 

DA 


dt 


)/A 


. . 2.8 
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Thus the eqn. 2.8 is an ordinary differential eqn. , which 
describes the evolution of time average value of U. 

The Fluxes H , H , H ,H are defined as 

AS BC CD BA 


H 


AB 


AB 


A Y 


AB 


AB 


A X 


AB 


H=FAy - GAX 

BC BC BC BC BC 


H=FAY - GAX 

CB CB CB CB CB 


H = F AY 

BA BA BA 


- G AX 

BA BA 


where F , F , G 

AB BC AB 


on sides AB , BC , CD , 


AX = X - X 

AB B A 


AY = Y - Y 

AB B A 


and so on respectively. 


, G are respectively values of F,G 

BC 

DA and where 


Let F. .be the mean value of F in elementary cell, see 


.J 


Fig C2.2). The straight forward way to define the fluxes is 


H = 4 C F. ^ . + 

AB 2 t-H ,j 


V ,J AB 2 ^ ,J ^ *3 AB 
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Uith analogous expressions for other fluxes. If uniform 
cartesian mesh is considered, these formulas yield the usual 
central differencing. In the same manner, it is possible to 
define a non-center ed scheme by assuming, for instance, that 
the mean of F on AB is defined as weighted average. 

* F + (1-61)* F instead of 4 ( F.^ .+ F. . ) 

j ; 2 L,j 


If u is positive then 6i = 1 
If u is negative then iii = 0 
for average value <.> =0.5 

Such a definition with 6i =1 or 61 =0.5 is used in Finite volume 


methods suggested by Mac-cormack and Paullay (1972). 
Ultimately eqn (2.8) reduces to 


dU 

dt 


+ H + H ) 

Ob bA 


/A 


So to get the flow field the above eqn is to be solved by 
numerical methods. In the present analysis we are using Runge 


Kutta method. 




Fig: 31. 
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Chapt,er- 3 

Method o£ solution 
The eqn to be solved is 


dU 

dt 


= C H + 

AB 


H 


BC 


+ H 


CB 


+ H > /A 

BA ^ ' 


For a particular geometry, a grid is generated over its 
whole domain. The values of all the variables are found at the 
center of each cell. Apply the unsteady method of choice to 
integrate the eqns forward in time subject to the inviscid flow 
boundary conditions of no flow normal to the surface and 
prescription of appropriate free stream conditions. In the 
present analysis we are using explicit shock capturing method. 

To solve the eqns we require the values of fluxes and 
the Areas of each cell, which are computed as follows. 


Me-Lhod bo find bhe area of IKe cell 

It is derived from the area under the curve method. 

Consider a cell as shown in fig( 3.1). Divide the cell into two 

triangles. Let -C * 2 *^ 2 ^ coordinates of the 

vertices of the triangle, then the area of the triangle (*, 2 , 3 ) 
given by 

A=0.5*Cx (y - y,)+x (y - y )+x(y - y )> ..(3.1) 

B=0.5*Cx^Cy3- ya» ..(3.2) 

Thus the area of cell is given by the sum of A + B .Thus by 

this the area of all the cells can be found. 
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Now the aecond step is to find the fluxes. 

pVL 

pVL^ + p 
puv 
puh 


li = pv + p 

pvh 

o 

Now let us represent them in terms of ( I,J,K),where ( I,J) are 
coordinate points and K=parameter. 

UC I.J.l) =p 
U (I,J,3) =puv 

Similarly the values of F and G can be written. The fourth 
parameter we are not considering as we have assumed constant 
enthalpy. So energy eqn does not come into account. Now to 
find u 

u = U UCI,J,1) 

Similarly v= U (I,J,3)/ U(I,J,1) 

So once p, u, and v are determined then p can be found by the 
eqn 


U = 


P 

pxi 

pv 

e 


F= 


and 


pv 

pVLV 

2 


P = 


P 

F 


[ R - ( y - 1 )( u* + V® )] 
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Thus by knowinfi the values of u ,v, p .and p the fluxes H can 
be calculated as defined earlier. Thus now applying proper 
boundary conditions the problem is solved. 

3.1 Test, case C Flow over a wedge 

PV3’. 3-2, 

The free stream conditions are 

1 . Mach number = 2 

2. Pressure = 1 bar 

3. r = 1.4 

4. R =287 J/ kg. K 

The unknowns to be found in the entire field 


1. 

P - density 


2. 

u - velocity 

components in x 

3. 

V - and 

y directions. 

4. 

p - pressure 



Uhen a supersonic flow passes over an compression corner 
a shock is formed. After the shock the Mach number decreases 
and the pressure, density and temperature rises 

discont inuously . The shock angle ft can be found from B-ft-U 
relations, Ref ( 4 ). Presently we have used Maccormack 
relationsto find the density, velocity, pressure, and Mach 
number after the shock. 

P=l, p =1, M>1 

Eqns are in non dimensional form.- 
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2y tt*. 3 i n fta h k ~ Cy ~ 1 ) 
__ 

— = (y^-DHo’^Sin^/Sahk 

*^0 Cj' + l)lto®Sin“^«hk + 2 



-l)Ifa Sin /?ahk +2 

”o Vsin C/3ahk -«*) 2j'lIo Sin/?ahk -(y-1) 

Computational mesh 
Dx = Ax = 1/16 

Dy = Ay = 2* Dx* (tan/? + u(|> + c<|> )/C u<j> + c<|> ) 
c<j> = J' Y .T 

u<J> = Moo * c<|) 

Thus the step length Ax and Ay are found. 

Now the time step is given by 

Dt= Smaller value of ( Dx , Dy ) 

u<(> +c<(> 

This is due to Cour ant , Friedrichs , and Levy (1928) 
restriction. 

BOUNDARY CONDITIONS 


.(3.1.1) 


.. (3.1.2) 


..(3.1.3) 


stability 


The boundary conditions involved in solving this problem are as 
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Sul^ersonic out^lo'jo 



pg: 3-B>- 
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follows : 

Inlet, exit, and solid wall boundary conditions are 
enforced during calculation of fluxes through the faces 
situated on the edge of the computation domain. For solid 
walls the mass flux through the wall vanishes ( q.n =0 ).The 
only contribution to the wall flux is the pressure acting on 
the wall, which is obtained by the second order extrapolation 
using the averages and the gradients of the cells adjacent to 
the walls. 

Fluxes through the inlet and the exit planes are 
obtained using ghost points outside the computational domain. 
For supersonic intakes the ghost point values are set at inlet 
conditions ,where-as for supersonic exit all quantities are 
extrapolated from the interior. For subsonic inf low, entropy , 

total enthalpy and flow angles are specified, while pressure is 
extrapolated to the inlet plane. C Using averages and the 
gradients of the first cell inside the domain). For subsonic 
outflow static pressure is specified and entropy, total, 
enthalpy and flow angle are extrapolated. Ref. (iO)- 
Mathema-bical representation of B .C 

There are four boundary conditions F‘3 C3*3^ 

1. Free stream boundary condition 

2. Uall boundary condition 

3. Two supersonic out flow conditions 
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1. Free stream B.C. 

For supersonic case fix all the variables i.e.Fix p, u, v, p 
For subsonic case fix Po ,To and every time calculate p, u, 
V, p . Flow at the inlet is of free stream case. Thus we know 
all the variables at the inlet and we can find the changes in 
flow as time step changes. 

2 . Uall B.C 

There is no flux in the direction normal to the wall. This 
is proved by an example. If we have a curved surface velocity 
in normal dirdirection is zero, i.e. V.n = 0 Let lx and ly be 
the direction cosines then ulx+ vly= 0. In our case the total 
flux is / F.nds where F= (F,G ) and F,G are defined 
earlier. So, 

Fix + Gly = pulx + pvly 
a pCulx + vly ) 

=0 as ( ulx + vly = 0 ) 

= ( pu*+ p )lx + C puv ) 

= pu(ulx + vly ) + plx 
= plx 

Similarly ( puv)lx +(pv*+ p)ly 

a ply 
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Thus we get 


F = 0 
F CI,J,2) = p 
F CI,J,3) = 0 


G = 0 
G CI,J,2) = 0 
G CI,J,3) = 0 


Thus at wall only pressure term is there and rest all the 
fluxes are zero. The pressure is found by linear 
extrapolation. 

Outflow B. C (at exit ) 

In case of supersonic flow it is just the linear extrapolation. 
From the central difference formula F. = F. + F. . 

V L — i I +1 

2 


F = 2* F. - F. 

t+l l L “i 

F = 2* F. - F. 

i i -i i -a 


Thus the value is extrapolated by using the values of two 

previous cells. 

outflow B.C (at top ) 

It is assumed that the outflow is supersonic and linear 
extrapolation is used. 


F. 

i 


F + F.^ 

J-4 J+i 


F. = 2* F. ^ - F. ^ 
j J-i 


Thus the fluxes at all points including the 


boundaries can be 


found. 
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RUMGE KUTTA METHOD 

After finding all the fluxes we substitute in the 
differential eqn. It being a simple first order differential 
eqn. Runge kutta method is used to solve it. Ue are using 
fourth order Runge kutta method for solving the differential 
eqn (3.1). It gives higher time step by y 2. times theoretically 
and also the error in computation reduces. The error graph 
falls sharply as the order increases , but above fourth order the 
changes are less compared to the computing time. 

The routine is written as: 

= u" 

U* = - 1/4 * At R* 

U® = U* - 1/3 * At R^ 

U® = - 1/2 * At R® 

U"* = U* - At R® 

U"'"* = U* 

Here and u'^^*are values at the beginning and the end ol 

t h . . . 

n time step. 
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CHAPTER 4 


DISSIPATION 

Aa we go for higher Hach numbers of the order of Mach 5 
it is found that there are many oscillations and solution 
instead of converging just blows up. In such applications it 
is often necessary to add an artificial viscosity (dissipation) 
term to the schemes described before. Such a term is 
necessary to; 

1. Select solution satisfying the entropy condition. 

2. Prevent nonlinear instabilities. 

3. To damp spurious oscillations of numerical solution. 

The concept of introducing an artificial term to smear 
out shocks was first introduced by von Neumann (1944). It was so 
devised that it was effective only in the regions where high 
compressions occur, i.e where shocks are building up. Away from 
the shock region ,the effect is negligible. The Euler eqns. are 
nondissipative hence it is desirable to have a scheme with 
minimum dissipation, for higher accuracy. 

In the present case there were many oscillations so some 
dissipation was added, to suppress the odd and even point 
oscillation. The damping added is proportional to difference 
between local value of enthalpy and the free stream value. 
Thus the Euler eqn becomes: 
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The dissipation operator which we used is a 
of second and fourth order differences. These 
separately in X and Y directions. 

Ex for density eqn 


Dp = D p + 

D 

X 

y 

D p = d. , 

- d 

u +iyz ,j 

i.-iy'a ,} 

D p = d. 

- d . 

y i,4+i^a 



combination 
are done 
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In the above eqns K is a constant .however the outcome 


is 
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Hl 

positive ,of the order ( x ), and proportional to 
difference of pressure. It does seem to do 
suppressing the oscillation around shocks. 


E. 

v j 


{ iZt 

= MAXC O.C K, - E. . )) 


i 

vhara K is second constant. 


the second 
a nice job 
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(::hapter 5 

RESUI.IS AND DISCUSSION 

WEDGE FLOW CASE 

Thl.H problem was taken as a trial case to test 

whether the algorittun is working well or not. This case is 
taken from Ref. (4) figure as shown in (3. 2). The initial 
condition is uniform flow. The wedge dimensions and the 
free stream conditions are described earlier. The wedge 
angle is 4,5i firstly a computational mesh is generated over 
the entire wedge and values of each variable are found at the 
node points, this goes on till the steady state is reached. 

The analysis was carried by two different grids, one, 
mesh along the shock wav^i and other perpendicular to 
the wedge, for Mach 2 flow.. According to the standard test 
data available when me.sh is aligned to the shock wave the 
density or the pressure variation along the wedge is as shown 
in fig (b B - b. 10) . ti l I. the shock the density is constant and 
sudden jump at the location of the shock and again constant 
over the entire length t>f the wedge. After computing the 
results, they were found satisfactory. 

When the grid was made perpendicular there were many 
oscillations in the shock region and, then uniform as the shock 
strength reduces. Usually this type of grid is used in solving 
the practica.L problems, whe.reas the former is the ideal case and 



useful for a particular problem where shock location and shock 
angle can easily be found out, i.e for simple geometries. 

For the wedge flow case the computation was done for 
two different cases: one w = 0.5 i.e. (central differencing) 
and other o> = 1 i.e. upwinding. In the former case there 
were wide oscillations (odd and even point) and the curve has 
saw tooth type profile. Fig(5.9). With upwinding the 
oscillations were damped to certain extent and results were 
slightly nearer to the ideal case. With upwinding the 
resulting finite difference eqn is only of first order 
accuracy , thus accuracy reduces slightly. 

To damp the oscillations further an artificial 
viscosity term (dissipation) was added along with upwinding 
and we got a smooth curve (ideal cu3rve).So due to upwinding and 
dissipation the oscil lations were fully damped .The results afe 
shown in graphs (5.6 - 5. 10). The curve trends were found to be 
same when tried for different Mach numbers. 

The intake geometry is as shown in the fig (1.1) and 
the size and the initial conditions are already specified in 
chapter 1. The physical domain was transformed into 
computational domain. A grid size of 55*10 was used in 
calculations. He grid was spaced equally and in same 
concent rati on everywhere 

To start with free stream conditions were assumed 
initially at all grid points except near boundaries wheire 



appropriate boundary conditions were used. Only the top 
boundary condition was changed i.e. wall boundary condition was 
introduced instead of supersonic outflow boundary condition 
compared to wedge case as discussed earlier, rest all boundaxry 
conditions remain same. 

The time step was determined by C.F.L restriction and 
integrated with time still steady state is reached. To access 
the quality of the solution the pressure profiles were 
compared at three locations on upper and lower walls. These 
lactations lie ahead, aft of, and in the separation region. 
Fig(i).4) and (5.5) show the pressuredistribution on upper and 
lower walls. When compared with Ajaykumar’s Ref.( 4.) solution 
the trend was found to be same but the numerical values were 
slightly differing. This may be mainly due to the efficient 
implicit and explicit ( combined )method used by him ,this 

method gives higher time step, secondly he used a 51 ♦ 51 

grid which gives higher resolution. Lastly he was dealing 
withaviscous problem so the amount of damping used may be 

different. He must have also used the grid refinement i.e 
close grids near the shocks. 

The pressure distributions at different stations were 
compared and the curve trends were found to be same. 

CONCLUSIONS 

With the above results we can conclude that that the 
results obtained by inviscld flow » 



th© inviscid flow solut-ioii first in order to establish the 
correctness of the algorithm which will form an integral part 
of general case as well as to have the preliminary results with 
much less effort in comparison with the general case. The 
extension of this analysis for the general realflow situation 
could not be taken up due to limited availability of time for 
thesis work. 
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